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ABSTRACT 

The temperature of the low-density intergalactic medium is set by the balance between 
adiabatic cooling resulting from the expansion of the universe, and photo-heating by 
the UV-background. We have analysed the Lya forest of eleven high-resolution quasar 
spectra using wavelets, and find strong evidence of a marked jump in the temperature 
at the mean density, Tq, of 60 ± 14 per cent over the redshift interval z = [3.5,3.1], 
which we attribute to reionization of Hen . The jump can be seen in all three of our 
spectra that straddle redshift 3.3, at a significance of > 99 per cent. Below z 3.1, our 
results are consistent with a smooth cooling down of the universe, as expected when 
adiabatic expansion dominates over photo-heating by a UV-background from QSOs 
and galaxies. We find no evidence of thermal fluctuations on scales > 5000 km s~^ 
larger than 50 per cent, which could be detected by our method, suggesting that the 
IGM follows a reasonably well-defined temperature-density relation. We demonstrate 
that the mean wavelet amplitude (A) cx I/Tq and calibrate the relation with hydro- 
dynamical simulations. We find Tq > 1.2 x lO^K at z > 3.6. Such high temperature 
suggest that Hi reionization occured relatively recent. 

Key words: cosmology: theory - intergalactic medium - hydrodynamics - large-scale 
structure of universe - quasars: absorption lines 



1 INTRODUCTION 

■ The intergalactic medium (IGM) can be observed in the 
' spectra of distant objects such as quasars through resonant 
absorption in the Lyman-a transition of neutral hydrogen 
(Bahcall & Salpeter 1965; Gunn & Peterson 1965). The 
availability of high resolution echelle spectrographs on large 
telescopes (HIRES on Keck & UVES on VLT) has provided 
us with data of unprecedented quality over recent years (see 
Ranch 1998 for a recent review). At the same time, a the- 
oretical paradigm within the context of the cold dark mat- 
ter cosmology has emerged (e.g. Bi, Boerner & Chu 1992), 

* The data presented herein were obtained at the W.M.Keck Ob- 
servatory, which is operated as a scientific partnership among the 
California Institute of Technology, the University of CaUfornia 
and the National Aeronautics and Space Administration. The Ob- 
servatory was made possible by the generous financial support of 
the W.M.Keck foundation. 

t Based on public data released from the VLT/UVES Commis- 
sioning and Science Verification and from the OPC porgramme 
65.0-296A (PI S.-D'Odorico) at the VLT/Kueyen telescope, ESO, 
Paranal, Chile. 



well tested by numerical hydrodynamical simulations(Cen et 
al. 1994; Zhang, Anninos & Norman 1995; Miralda-Escude 
et al. 1996; Hernquist et al. 1996; Wadsley & Bond 1996; 
Zhang et al. 1997; Theuns et al. 1998; Machacek et al 2000) , 
in which the absorption is produced by volume filling, photo- 
ionised gas, that contains most of the baryons at redshifts 
2 ~ 3 (see e.g. Efstathiou, Schaye & Theuns 2000 for a recent 
review). The absorbers are locally over dense extended struc- 
tures, close to local hydrostatic equilibrium (Schaye 2001). 
The combination of a predictive theory and superb data has 
led to a veritable revolution in IGM studies. 

On large scales, the gas distribution is similar to that 
of the underlying dark matter, but on small scales pressure 
forces smooth the gas distribution, erasing most of the small 
scale power. Combined with thermal broadening and pecu- 
liar velocities, this results in a strong suppression in the 
amplitude of the flux power spectrum on scales < 200 km 
s^^ (e.g. Theuns, Schaye & Haehnelt 2000). Another conse- 
quence of the thermal content of the gas is a cut-off in the 
distribution of line widths 6 below ~ 20 km s~^ (Schaye et 
al 1999) , which is clearly present in observed line samples 
as well (e.g. Kirkman & Tytler 1997). Both these signatures 
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can be used to infer the IGM temperature by calibrating the 
cut-off using simulations (Schaye et al 1999, 2000; Ricotti, 
Gnedin & ShuU 2000; Bryan & Machacek 2000; McDonald 
et al 2000). 

The IGM temperature can be used to constrain the 
reionization history because the thermal time scales are 
long in the low density gas probed by the Lya absorption 
(Miralda-Escude & Rees 1994). Consequently the gas re- 
tains a memory of its heating history. At over densities 
Sp/{p) < 3, shock heating and radiative cooling are unim- 
portant and photo-heating of the expanding gas introduces a 
tight density-temperature relation T/To — {p/{p)y'~^ (e.g. 
Hui and Gnedin 1997). A sudden epoch of reionization will 
make the gas nearly isothermal (7 ~ 1) at a temperature of 
~ lO^K. Away from reionization, the density-temperature 
relation steepens again, asymptotically reaching 7 ~ 1.6. 
The steepening occurs because denser gas has a higher neu- 
tral fraction, and hence a larger photo-heating rate. If the 
heating rate can be computed reliably, given the observed 
sources of ionising photons, then a measurement of the en- 
tropy of the gas can in principle determine the reioniza- 
tion epoch. If the sources responsible for reionizing hydrogen 
have a soft spectrum, there may be more than one reioniza- 
tion epoch since helium reionization requires harder pho- 
tons. 

There is some observational evidence that Hell 
reionization occurs around 2: ~ 3, based on the occurrence 
of large fluctuations in the Hen optical depth (Davidsen et 
al. 1996; Reimers et al. 1997; Heap et al. 2000), and the 
change in hardness of the ionising spectrum as inferred from 
metal-line ratios (Songaila & Cowie 1996; Songaila 1998). 
The interpretation of the optical depth data is not straight- 
forward however (Miralda-Escude et al. 2000), and there is 
some discussion in the literature about the reality and in- 
terpretation of the claimed change in metal-line ratios (Bok- 
senberg, Sargent & Ranch 1998; Giroux & ShuU 1997). Note 
that, if Hen reionization does not occur instantaneously, var- 
ious measures of the Hen abundance change may well find a 
range of 'reionization epochs', given that they may sample 
different over densities. 

Reionization is not expected to be instantaneous be- 
cause ionisation fronts expand much faster into low density 
voids than into the higher density filaments (e.g. Gnedin 
2000a). If the sources of ionising photons are bright but 
scarce, then reionization will not be complete until the large 
ionised bubbles overlap. A consequence of such patchy reion- 
ization is that different regions might follow different p — T 
relations and so have different values of Tq. Such 'temper- 
ature fluctuations' may lead to apparent discrepancies be- 
tween the different algorithms discussed earlier to measure 
To and 7, since different methods may give other weights to 
regions in the same spectrum. 

The aim of this paper is to test whether the data are 
consistent with a single well defined p — T relation and to 
investigate whether Hen reionization can be detected by a 
jump in temperature. We use the wavelet analysis proposed 
by Theuns & Zaroubi (2000), and apply it to a set of high 
resolution QSO spectra spanning a wide range of redshifts. 
We use hydrodynamical simulations to demonstrate that our 
method can detect variations in the amplitude To of the 
temperature-density relation p/{p) oc (T/To)^~^ of order 50 
per cent, even when these occur on scales as small as ~ 5000 



km s~^ (1000km s"'^ corresponds to ~ 9.3 co-moving Mpc 
h~^ and a redshift extent Sz — 1.33 x 10~^ at redshift z = 3, 
in the currently popular flat, cosmological constant domi- 
nated model with matter density Qm = 0.3.) The method 
presented here does not use simulations to identify regions 
of different temperatures. Recently, Zaldarriaga (2001) ap- 
plied a similar analysis to QSO 1422-1-231 and showed that a 
model with two temperatures, that occupy comparable frac- 
tions of the spectrum, is constrained to have temperature 
variations smaller than a factor of 2.5. 

This paper is organised as follows. Section 2 presents 
the data and Section 3 the simulations used to demonstrate 
the method. Section 4 describes the wavelet decomposition 
and the statistical tools making use of hydrodynamical sim- 
ulations to illustrate the approach. Section 5 presents the 
results, and these are discussed in Section 6. Section 7 sum- 
marises. Readers not interested in the details of the method 
can skip directly to section 4.3.1, which provides an illustra- 
tion of the method for a simulated spectrum with imposed 
temperature fluctuations. 



2 THE DATA 

In our analysis, we have used data obtained with the high 
resolution spectrograph on the Keck I telescope (Vogt et al. 
1994), and UVES, the ultraviolet echelle spectrograph on the 
VLT (Kueyen) telescope (D' Odorico et al. 2000). We have 
combined line lists taken from the literature with publically 
available spectra (Table|), to obtain a set of eleven high 
resolution QSO spectra that span a wide range of emission 
redshifts Zcm = 1.7 — » 3.7. High resolution data is required 
since the lines are intrinsically narrow ~ 20 km s"^ The 
standard data reduction is described in the original refer- 
ences, with the exception of APM 0827-1-5255 which was re- 
analyzed after calibration of the flux scale to remove echelle 
order mismatches (Tzanavaris & Carswell 2002). All spectra 
have signal to noise of 40-50 per resolution element, and a 
similar resolution, R ~ 40, 000. Only two out of eleven of 
these spectra (QSOs APM 08274-5255 and 1422-1-293) were 
used in the analysis of the thermal evolution by Schaye et 
al (2000). 

Line profile fitting has been performed on the data. 
Voigt profiles are fitted to the absorption lines using a 
minimisation procedure, producing a list of lines with given 
column density, A'jjj (cm~^), width b (km s~^) and ab- 
sorption redshift z. Where possible, metal lines have been 
identified from line coincidences. This is in fact a crucial 
step, since such lines tend to be narrow and so might be 
mistaken for a cold Lya cloud. 

For QSOs Q0636-F680, Q0956-H122 and Q0014+813, we 
use the published line list given by Hu et al (1995), who used 
their own automated line fitting programme. All other spec- 
tra have been analysed with the semi-automatic line-fitting 
programme VPFn]T| (Webb 1987; CarsweU et al. 1987). The 
mock spectra from our simulations, described in the next 
section, are also fitted using VPFIT. 

The line parameters from Voigt profile fitting are not 

t for rintail.s on VPFTT. see 

http: / /www. ast.cam.ac.uk/~rfc/vpfit. htm] 
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unique (Kirkman & Tytler 1997). However, in the analy- 
sis described below, we never use the detailed properties of 
the lines, but rather analyse spectra as reconstructed from 
the line list. For QSO HE 1122-1648, we have redone our 
analysis using another line fitting programme, FITLYMAN 
(Fontana & Ballester 1995), and the results are nearly iden- 
tical. 



3 HYDRODYNAMICAL SIMULATIONS AND 
MOCK SPECTRA 

We use hydrodynamical simulations to illustrate the method 
described below. Briefly, these are simulations of a flat, 
vacuum-energy dominated cold dark matter model (matter 
density Qrn ~ 0.3, baryon fraction flth^ = 0.019, Hubble 
constant Ho = lOO/i km s~^ Mpc^^, h — 0.65 and nor- 
malisation (78 = 0.9. We have used CMBFAST (Seljak & 
Zaldarriaga 1996) to compute the appropriate linear transfer 
function. 

The simulation code is based on HYDRA (Couchman, 
Thomas & Pearce 1995) and combines Smoothed Particle 
Hydrodynamics (SPH, Gingold & Monaghan 1977; Lucy 
1977) with adaptive P3M gravity (Couchman 1991). It has 
been modified extensively by one of us (TT) in order to be 
able to simulate the IGM. The code has been tested compre- 
hensively, and has been parallellized in the OpenMP stan- 
dard. Non-equilibrium radiative processes (cooling, photo- 
ionization heating by a UV-background and Compton cool- 
ing off the cosmic micro-wave background radiation) are in- 
cluded with rates given as in Theuns et al. (1998). The evo- 
lution of the UV-background is chosen to give the thermal 
evolution as determined by Schaye et al. (2000, the simula- 
tion is referred to there as the 'designer simulation'). In this 
simulation. Hi and Hel reionize at a redshift 2 ~ 7 by soft 
sources, whereas Hell reionization is delayed to a redshift 
z ~ 3.4. The simulation box is 12.5/i~^ co-moving Mpc on a 
side (corresponding to ~ 1400 km s~^ at redshift z = 3), and 
gas and dark matter are represented with 2 x 256'^ particles 
of masses 1.45 and 8.25 x IO^Mq respectively. The resolu- 
tion of this simulation is sufficient to resolve the various line 
broadening mechanisms (Theuns et al. 1998), which is of 
course crucial for this type of analysis. At the same time, ef- 
fects of missing large scale power on the statistics analysed 
here, are not very important at sufficiently high redshifts 
> 2 (Theuns & Zaroubi 2000). 

In order to study the effect of temperature fluctua- 
tions on the properties of the absorption lines produced, 
we impose specific temperature-density relations, T/To = 
(p/(p))^~^, on our simulation outputs. We impose such 
a T ~ p relation on all gas particles with p < (p) or 
T(p) < 2 X To{p/{p))'<-^ and p < 10{p). Models C and 
H (for Cold and Hot respectively) have To = 1.5 x 10*K and 
2.2 X IO^'K respectively, with the same slope, 7 = 5/3. We 
have used these two models to make mock spectra of length 
~ 50000 km s~^, typical of observed Lyo? spectra of redshift 
~ 3 QSOs. (These combined spectra use sets of Voigt profiles 
pertaining to fits of spectra of individual sightlines through 
the simulation box. Each spectrum is fitted such that the 
maximum transmission occurs at the edges, by taking ad- 
vantage of the fact that the individual chuncks are periodic.) 
In order to make the comparison to data more realistic, we 



Table 2. Mock spectra used. (To,7) refers to the imposed 
temperature-density relation, y = 5 X 10^ km s~^ is the length 
of each spectrum. 



Spectrum 




SI 


(1.5 X 10'',5/3) for V < 0.5V 




(2.2 X lO*, 5/3) for v > 0.5V 


S2 


(2.2 X 104,5/3) for |i; - V/2\ > 2500 




(1.5 X 10*, 5/3) else 


S3 


(2.2 X 10", 5/3)1 



1 Mean flux F = 0.64 (0.6) for v < V/2 {v > V/2). 



impose the observed mean absorption on the spectra. In ad- 
dition, we also add noise and 'instrumental broadening' fol- 
lowing the procedure described in detail in Theuns, Schaye 
& Haehneh (2000). 

Below we analyse several such mock spectra (Table 
Spectrum SI is composed of spectra drawn from model C 
for its first half, and from model H for its second half. Spec- 
trum S2 is drawn from model H, except for a cold gap of 
length 5000 km s"^ drawn from model C in the middle of 
the spectrum. Finally, S3 is a single temperature model (H) , 
but the mean absorption is scaled to be respectively 0.64 and 
0.6 in the lower and upper redshift halves of the spectrum 
respectively. 

To illustrate the effect of varying 7, we also compute 
spectra for models with imposed T — p relations of (To, 7) = 
(1.5 X 10*, 1) and (2.2 x 10*, 1). (The neglect of extra broad- 
ening due to pressure effects and peculair velocities in spec- 
tra with imposed equation of state artificially decreases the 
differences in line shapes for different To models.) In addi- 
tion to these simulations, we use another set of four simula- 
tions evolved with 2 x 64'' particles in a box of size 2.5h~^ 
Mpc and the same ACDM cosmology as the others. These 
simulations have been run with different heating rates, to 
give them different values of logioTo = [4.0, 4.1, 4.2, 4.3] but 
the same value of 7 « 1.4. The mock spectra from these 
simulations are computed in the same way as the others. 

These mock spectra are then scaled to have the same 
mean absorption as the data. We have used VPFIT to pro- 
duce line lists for all mock spectra, and this allows us to treat 
data and simulation in an identical fashion for the wavelet 
analysis presented in the next section. 



4 METHOD 

4.1 Wavelet projection 

A discrete wavelet is a localised function with a finite band- 
width (see e.g. Press et al. 1992 for an introduction and 
original references to the application of wavelets to a wide 
variety of problems). This makes wavelets useful for char- 
acterising line widths in a spectrum, since the amplitude of 
the wavelet will be related to the width of the line, and the 
position of the wavelet to the position of the line. Pando 
& Fang (1996) used a wavelet analysis to describe cluster- 
ing of Lya absorption lines. Theuns & Zaroubi (2000) used 
the Daubechies 20 wavelet (Daubechies 1988) to characterise 
temperature fluctuations in a Lya spectrum. Independently, 



© 0000 RAS, MNRAS 000, 000-000 



4 Theuns et al. 

Table 1. QSO spectra used 



\alO\J 






AA 


^Lya 


C o HI III G nt s 


Ref 




1.719 


14.9 


3080- 


-3270 


1.53- 


-1.69 


VT,T /TTVFS 


\ 


J2233-606 


2.238 


17.5 


3400- 


-3890 


1.80- 


-2.20 


VLT/UVES 


1,2 


HE1122-1648 


2.400 


17.7 


3500- 


-4091 


1.88- 


-2.37 


VLT/UVES 


2 


HE2217-2818 


2.413 


16.0 


3510- 


-4100 


1.89- 


-2.37 


VLT/UVES 


1,2 


Q0636+680 


3.174 


16.5 


4300- 


-4900 


2.54- 


-3.03 


Keck/HIRES 


3 


Q0302-003 


3.286 


18.4 


4808- 


-5150 


2.96- 


-3.24 


VLT/UVES 


2 


Q0956+122 


3.301 


17.8 


4400- 


-5000 


2.62- 


-3.11 


Keck/HIRES 


3 


Q0014+813 


3.384 


16.5 


4500- 


-5100 


2.70- 


-3.20 


Keck/HIRES 


3 


Q1422+231 


3.620 


16.5 


3645- 


-7306 


2.91- 


3.60 


Keck/HIRES 


4 


Q0055-269 


3.655 


17.9 


4852- 


-5598 


2.99- 


-3.60 


VLT/UVES 


2 


APM08279+5255 


3.911 


15.2 


4400- 


-9250 


3.20- 


-3.72 


Keck/HIRES 


5 



B-band magnitudes from the SIMBAD astronomical database, except for QSO APM082794-5255, for which the i?-band magnitude is 
from Irwin et al. (1998). 

References: 1. Kim et al. (2001a); 2. Kim et al. (2001b); 3. Hu et al. (1995); 4. Rauch et al. (1997); 5. Ellison et al. (1999) 



Meiksin (2000) used the same wavelet basis to show that 
most of the information in a spectrum is carried by a small 
fraction of the wavelets. More recently, Jamkhedar, Bi and 
Fang (2001) used a multi-scale analysis based on wavelets 
to better describe the transmission power spectrum taking 
into account uncertainties in the applied continuum. 

Discrete wavelets are also a set of orthogonal basis func- 
tions. Consequently, we can borrow notation from quantum 
mechanics to write any function as a sum over its wavelet 
projections, \ip) = ^ {i/>n|i/') \'<Pn)- Here, is the wavelet 
function of level n, {ipn\ip) the projection of \^p} onto \tpn}, 
and n the 'quantum number' of the basis function. Charac- 
teristic for wavelets is that they have a finite extent in both 
the spatial direction, and the frequency domain, i.e. they are 
like a wave-packet instead of a single momentum wave. So 
we'll need two indices to denote a given wavelet basis func- 
tion, lipn^p), where n denotes the characteristic frequency 
('width'), and p the position of the wavelet. 

In this paper we will investigate the properties of the 
spectrum in terms of a given single frequency of the wavelet, 
i.e. we will study 

IV-W) = ^(^n,piV)l^n,p). (1) 

P 

Note that summing |i/'(n)) over n will give back the original 
function \tp}. We are left to tune the frequency n to make it 
most sensitive to the temperature. We have chosen to use the 
Daubechies 20 wavelet, as in Theuns & Zaroubi (2000), who 
give examples of the decomposition of a spectrum in terms 
of that wavelet. Note that computing the wavelet projection 
Eq. (jl|) is computationally very similar to performing a fast 
Fourier transform, and many standard computer packages 
have wavelet transforms built-in. The usage of wavelets is 
not very common in astronomy, we would therefore like to 
note that the computation of the projection Eq. (0) requires 
a handful of programming lines in standard software pack- 
ages such as IDL, and much less than a second of computer 
time to evaluate. For a given choice of the wavelet, the de- 
composition Eq. (|l]) is unique. 

Theuns & Zaroubi (2000) showed that the projection in 
Eq. (Fl) is already very sensitive to the temperature of the 



gas, because narrow lines tend to generate larger wavelet 
coefficients than broader lines, for a suitable choice of n. 
Here we describe four improvements to that method. 

(1) Since the \tpn,p) are orthogonal basis functions, they 
change sign at least once, in contrast to an absorption spec- 
trum. In particular, the Daubechies 20 wavelet does not 
look like an absorption line at all, but is more similar to 
its derivative. So we found it advantageous to project the 
scaled derivative of the spectrum F{\) with respect to ve- 
locity V, dF/dv{F + rj)~^ . Here, the 'velocity' v is defined as 
a function of wavelength A through dX/X = dv/c, where c is 
the speed of light, and the parameter rj = 0.2 is introduced 
to avoid division by zero close to the zero level. 

(2) The positional localisation of the narrow wavelet (width 
~ 15 km s~^) that we will use below is far better than what 
we require to determine a jump in temperature. We take 
advantage of this to improve the frequency resolution - and 
hence the sensitivity to small changes in the line widths - of 
the wavelet projection at the expense of its localisation by 
computing 

\ip{n,5)) = y^^rnax(^_i/2fi,i/2)xs{\{i^{n,p)\ip{~5)) 

V 

X \^(n,p + 5))\)- (2) 

So we shift the input function by 5 to the left, |7/)(— 5)), 
project it onto the wavelet, and shift the result by the same 
5 to the right. We do this for 5 = -1/2,0,1/2 times the 
wavelet's width, and at each position take the maximum 
absolute value of the three projections, and denote the result 
by \^{n, 5)). This step is not crucial to the method, but does 
improve its sensitivity in finding temperature fiuctuations 
imposed in the mock spectra. 

(3) Observed spectra contain metal lines in addition to 
hya lines. As described in the data section, we use VPFIT 
to fit Voigt profiles, V{bi, Ni, Xi), to the spectrum. Here, 
V{bi, Ni, Xi, Xi) denotes a Voigt profile of species Xi with 
column density A'^;, width bi and centred on Xi. Using this 
decomposition, we can reconstruct an (almost) 'metal free' 
spectrum. 
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FvP = exp{-Y^Vib„N„\„Lya)). (3) 



(Note that we still impose instrumental broadening on this 
reconstruct ion.) 

The decomposition Eq. (0) is very good at identifying 
very narrow lines < 15 km s~ . When applying the method 
to data, it is sometimes impossible to decide whether such a 
narrow line is a metal line or not, because the corresponding 
other metal transitions may not fall in the observed region 
of the spectrum, or fall on top of a strong Lya line. In cases 
where a strong wavelet signal results from a single Voigt pro- 
file, we have decided to flag the line as a potential metal line, 
and not consider it any further in the rest of the analysis. 
Typically only a hand full of lines (out of several hundreds) 
are removed per spectrum. This is of course a somewhat sub- 
jective procedure, but is unlikely to introduce any positive 
detection of a spurious correlation signal. 

Using Fyp instead of the original spectrum F has the 
added beneflt that it is also noise free, making it easier to 
take the derivative (step 1), and will also be extremely useful 
in the statistical analysis. The decomposition in Voigt pro- 
files is not unique, but since we only use the reconstructed 
spectrum, and not the detailed properties of the fitted lines, 
we expect the spectrum Fyp to be nearly independent of the 
actual manner in which Voigt profiles were fitted. In partic- 
ular, we were unable to detect any systematic differences in 
line shapes between the original spectrum and the recon- 
struction, using either our mock spectra or observed data. 
We also redid the analysis for one QSO spectrum (QSO 
HE 1122-1648) using FITLYMAN (Fontana & Ballester 
1995), another line fitting programme, and found the same 
results. 

In the absence of evolution in the b parameters, line 
widths AA will increase with redshift oc {1 + z). This pure 
'expansion' effect can be compensated for by rebinning the 
spectrum to velocity space, and VPFIT performs such a 
rebinning. 

(4) The spectrum also contains strong lines. Their square 
shape tends to generate large wavelet amplitudes, an effect 
similar to the Gibbs phenomenon familiar from Fourier anal- 
ysis. In addition, in the saturated region of the strong lines, 
the wavelet amplitude becomes zero. These fluctuations in 
the wavelet amplitudes do not contain any information on 
the IGM temperature and so we ignore them. We begin by 
identifying all regions of the spectrum within 1 per cent of 
being black. All pixels within 1.5 times the wavelet-width of 
such a region are ignored in the statistical analysis described 
in the next section. 

In summary: we start by fitting the original spectrum 
F{\) with Voigt profiles, and reconstruct the spectrum from 
the fit, using only those lines not identified as (or sus- 
pected to be) metal lines. For this spectrum, Fvp{v), we 
compute the wavelet projection of its normalised derivative, 
Fvp{S) = J2p'^"'^(-i-/2,o,i/2)x5{\{tp{n,p)\dFYp/dvx{Fyp + 
0.2)~^{-S))\^p{n,p + 5))\). We use the correlations in Fyp{5) 
to quantify temperature variations along the spectrum. Fig- 
ure y illustrates the procedure. It shows two stretches taken 
from spectrum SI of length 5000 km s~^, corresponding to 
model C and the 50 per cent hotter model H, respectively. 
The wavelet amplitudes tend to be significantly larger in the 
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Figure 2. Cumulative distribution of wavelet amplitudes for 
spectrum SI (full line labelled all), its hot half (long dashed line) 
and its cold half (dotted line). Vertical lines indicate the maxi- 
mum differences between the cumulative distributions. 




Figure 3. Maximum difference A between the cumulative distri- 
bution of j4 in a window of size 5000 km s~^ starting at velocity v, 
and the cumulative distribution of A over the whole spectrum, for 
the simulated spectrum SI. Negative values of A indicate regions 
where the window has a larger fraction of high wavelet amplitudes 
than the spectrum as a whole. A ~ —0.2 for the first half of the 
spectrum, and ~ -|-0.2b for the second half, showing that A is 
able to recognise that the first half of spectrum SI is cold, and 
the second half is hot. 

spectrum taken from the colder model. The wavelet ampli- 
tudes can therefore serve as a thermometer. 

4.2 Statistical analysis 

Jf..2.1 Cumulative distribution of wavelet coefficients 

In the previous section we illustrated the strong depen- 
dence of the wavelet amplitude on the temperature of the 
gas — cooler gas generates a larger fraction of narrow 
lines with correspondingly higher wavelet amplitudes A, 
on average. So we can now recognise temperature fluctu- 
ations by identifying regions over which the wavelet ampli- 
tudes are unusually large (a cold region) , or unusually small 
(a hot region). We will quantify the extent to which the 
wavelet amplitudes in two regions of the spectrum differ, 
by A{vi,dvi;v2,dv2), the maximum difference between the 
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Figure 1. Flux (lines, left hand scale) and wavelet amplitude averaged over 100 km s-l (filled histo gram, right hand scale) for a stretch 
of spectrum of length 5000 km s~^ from simulation C (To = 1.5 X lO^K, top panel) and the hotter simulation H (To = 2.2 X lO^K, 
bottom panel). The wavelet amplitudes tend to be significantly larger in the colder model. 
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Figure 4. Stretch of spectrum of QSO 1422+231. A piece of spec- 
trum between two crosses is identified as an 'absorption feature' 
in the randomisation procedure. Voigt profile fitted absorption 
lines with central wavelengths that fall within the same absorp- 
tion feature always keep their relative spacing when randomising 
the spectrum, thereby keeping the detailed line shape identical 
between original and randomised spectrum. 
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two cumulative distributions of A in the regions [iii, 
and [v2,V2 + dv2\ (velocities v are in km s""*^). (Note that 
/^{v-i_,dvi\v2,dv2) = — A(i!2, 9v2; fi, cJ'Ui)-) We will denote 
this difference as /S.{v,dv) when one of the windows refers 
to the whole spectrum. 

The cumulative distribution P(< A) of the wavelet 



Figure 5. Histogram of the lengths of absorption features iden- 
tified in QSO 1422+231. 



amplitude A for the two windows shown earlier in Fig. ^ 
([1.9 x 10^,5 X 10=*] and [3.55 x 10"*, 5 x 10^] respectively) is 
plotted in Fig. ^, together with the cumulative distribution 
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Figure 6. Probability distribution P(A) for windows of size 5000 
km s~^ for windows drawn from 200 randomised spectra of spec- 
trum SI (dotted line). Dashed line refers to P(A) for the original 
spectrum SI, but only includes window pairs that fall in the same 
half of SI, thick full line refers to window pairs that fall in differ- 
ent halves. The dotted and dashed histograms are very similar, 
because both refer to single temperature models. The histogram 
referring to different temperature window pairs (full line) differs 
significantly from the case where both windows refer to the same 
temperature. This type of statistic is therefore useful in recognis- 
ing regions with different temperatures. 



of A over the whole spectrum of SI. The three distributions 
appear to be quite different, as expected. The differences A 
between the cumulative distributions are shown as vertical 
lines. Note that Theuns & Zaroubi (2000) also used P{< A) 
to characterise the T — p relation. 

The quantity A(u, dv) can be used as a (uncalibrated 
as of yet) thermometer. In Fig. ^ we plot A(w,5 x 10'') as 
a function of the starting position v of the window for the 
simulated spectrum SI. Notice that as long as v is in the first 
half of the spectrum - in which case the window is drawn 
from the cold simulation C - A ~ —0.2. There is a sudden 
transition around v ~ 2.5 x lO** km s~^ where A jumps to 
values ^ -1-0.2, as the window starts to fall in the stretch 
of spectrum drawn from the hotter model H. We want to 
stress again that, for a given wavelet basis, A(«, 5 x 10^) 
follows uniquely from projecting the spectrum on a set of 
basis functions. 

Our aim is now to use A (ui , dvi ; V2 , dv2 ) to find regions 
with different temperatures, and use A{v, dv) as a measure 
of To. However, how can one judge whether a given value of 
A is statistically significant? The usual way (Kolmogorov- 
Smirnov test) to decide whether two data sets are drawn 
from the same underlying distribution is to evaluate a par- 
ticular function of both [A], and of the effective number of 
degrees of freedom A''. Unfortunately, in the present case it 
is unclear how to determine A^, since pixels are correlated. 
One way to make progress is to use randomised spectra, in 
which any temperature fluctuations have been destroyed in 
the randomisation procedure, to calibrate the distribution of 
A and hence determine how statistically significant the dif- 
ference seen in Fig. ^ is. Note that, if a randomised spectrum, 
by chance, produces a a peculiar value of A, it is not due to 
temperature fluctuations, since the temperture-density rela- 
tion is not different from that of the rest of the spectrum. 
The full procedure to assign a statistical significance to un- 
usual regions is described next. 



4-2.2 Randomised spectra 

The aim of this procedure is to produce new spectra from 
the data, in which the absorption lines have the same shapes, 
but any correlation between the lines is destroyed. One could 
in principle randomise the positions of the Voigt profiles 
in Eq. (^, but the resulting spectra turn out to be quite 
different from the original one. The reason is that many 
absorption features are composed of several Voigt profiles 
(often with large error estimates for the fitted parameters), 
some of which may be quite narrow and are introduced to 
obtain a good fit. When randomising Voigt profiles, such a 
narrow line will tend to occur on its own, and this causes 
clear systematic differences between the randomised spectra 
and the original one. 

We have chosen instead to randomise the positions of 
the absorption features ('absorption lines' as opposed to 
Voigt profiles, in what follows) themselves. An absorption 
line is defined as a stretch of spectrum between two local 
maxima in Fvp. Voigt profiles are now assigned to a line, 
if the centre of the profile falls between the two maxima 
that define the line. Randomised spectra are then obtained 
from randomising the positions of the lines, making sure that 
lines do not overlap (but the Voigt profiles may) , resampling 
them with replacement. In the tests presented below, we will 
show that these randomised spectra are not distinguishable 
from the original one, in the absence of imposed temperature 
fluctuations. 

An example of the identification of lines for part of the 
spectrum of QSO 1422-1-231 is shown in Fig. ^, where local 
maxima are indicated by crosses. We use only those local 
maxima where the flux F > 0.5. A histogram of the widths 
of lines is shown for reference in Fig. ^. The typical widths 
of the lines is ~ 100 km s~^, with a tail to larger values. 

Given these randomised spectra, we can compute the 
probability distribution -P(A), by randomly sampling ~ 
2 x lO'* window pairs for each spectrum. The result is plot- 
ted in Fig. P (dotted line). Superposed is -P(A) for the orig- 
inal spectrum SI, for window pairs drawn from the same 
half of the spectrum (dashed line) and from different halves 
(fuU line). The randomised spectra have a similar -P(A) as 
the original spectrum SI, when both windows are drawn 
from a single temperature region, but the mixed tempera- 
ture -P(A) (full line) is very different. The probability dis- 
tribution -P(A) from the randomised spectra allows us to 
assign a statistical significance to a given value of A. 

4-2.3 Cluster analysis 

In the previous sections we demonstrated how the difference 
A can be used as a thermometer, and how the statistics of 
P(A) can be used to decide whether a difference in temper- 
ature between two windows is statistically significant. Here 
we will combine these two results to identify regions in the 
spectrum that have a statistically significant different tem- 
perature as compared to the rest of the spectrum. 

In order to do so, we study correlations in the fleld 

Q{vi,V2;dv) — A{vi,dv)\A{vi,dv;v2,dv)\ . (4) 

Recall that large values of \A{vi,dv;v2,dv)\ suggest differ- 
ent temperatures between the two windows [vi , vi -\- dv] ver- 
sus [v2 , V2 + dv] - Because we multiply this quantity with 
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Figure 7. Contour plot for the correlation field Q for spectrum 
SI. Contours are drawn at levels 0.02, 0.05 and 0.09 for positive 
values (full lines) and negative values (dashed lines) . Window size 
dv = 5000 kms^^. Contour level -0.021 delineates the cold region, 
and level 0.021 the hot region. 
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Figure 8. Same as Fig. ^ for spectrum S2. The cold gap in the 
middle of the spectrum stands out. 



/S.{vi,dv) to obtain Q, we can also judge whether window 
[vi , vi + dv] is cold or hot (in which case Q ^ and Q 3> 0, 
respectively), as we demonstrated in Fig 

A contour representation of Q is plotted in Fig. ^ for 
spectrum SI, which shows that Q < in the top left hand 
corner, and Q > in the bottom right hand corner. Recall- 
ing our previous results, this indicates that the first half of 
the spectrum is cold, and has a significantly different tem- 
perature from the second half, which consequently is hot. 



Recall that this is indeed the case for our simulated spec- 
trum SI. 

The next step is to determine the extent of the region 
where the temperature differs significantly from the rest of 
the spectrum. We do this using a cluster analysis as fol- 
lows. Given the normalised probability distribution P{A), 
determined from randomised spectra without fluctuations, 
we determine A/ such that J^^ P(A)dA = /, for / = 0.1 
and / = 0.9. That is, only 10 per cent of window pairs 
have A < Ao.i, and equally ten per cent of pairs have 
A > A0.9. We define Qi = |Ao.i A0.9I as a good first in- 
dicator of whether a window has a peculiar temperature, by 
comparing Q with Qi. Windows with Q < —Qi are likely 
to be peculiarly cold, those with Q > Qi peculiarly hot. In 
Fig. M, the lowest contour level shown is Qi, and indeed, 
these contour levels delineate the low (cold) and high (hot) 
regions of the spectrum. 

Another example of the field Q, now for the spectrum 
S2, is shown in Fig. ]A Recall that S2 has a single cold region 
located in [2.2 x 10^,2.7 x 10**]. This pecuhar region falls 
below the contour level — Qi and is therefore neatly detected 
by the procedure described so far. Note that there are other, 
smaller, regions where \ Q\ > Qi in this example. 

We now define a cluster as a connected region where 
\Q\ > Q;. In practise, we identify clusters by interpolating Q 
onto a grid of grid size dv/5 using cloud-in-cell interpolation. 
A cluster is then a set of connected grid cell^, which each 
have \Q\ > Qi. The weight of a cluster is the integral of Q 
over its cells. The motivation for doing this, is that a region 
of unusual temperature will have a large extent in a plot 
such as Fig. ^, and so will have a large weight. We will use 
the randomised spectra to judge the statistical significantly 
of those weights. 

When identifying clusters, we impose 'periodic bound- 
ary conditions' in the vertical direction (i.e., cells at the up- 
per and lower extreme, but the same horizontal position vi , 
are also neighbouring cells). This may seem strange at first, 
but consider the cluster identified in Fig. ^. Without periodic 
boundary conditions, there would be two clusters (of roughly 
half the size) above the cold region [2.2 x lO'', 2.7 x lO'']. The 
reason that these two clusters are not connected, but appear 
to have a hole around V2 ~ 2.5 x 10** is of course because 
there, v\ ^ V2, and both windows refer to the same stretch of 
spectrum - hence Q ~ 0. So clearly that region of the plane 
is not peculiar as it does not refer to windows of different 
temperatures. However, if the cold region had been closer to 
the start, or the end of the spectrum, there wouldn't have 
been a hole, and hence the cluster had been identified with 
nearly twice the weight it has now. This is clearly not what 
we want. Imposing periodic boundary conditions connects 
these two clusters again into one. In this case, the weight of 
the cluster does not depend on its position. Note that clus- 
ters that correspond to low temperatures have a negative 
weight. 

After cluster identification, we can introduce the new 
field, C, defined on the grid, where C{vi,V2) = if the cell 
does not fall into a cluster, and equals the weight of the 
cluster, if it does. We finally project C onto the spectrum. 



3 Grid cells that neighbour each other either horizontally or ver- 
tically, are called connected. 
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Figure 10. Temperature T1.2 = Tq X 1.2''"^ at an overdensity 
1.2, versus the inverse {A)~^ of the wavelet amplitude for models 
with different values of Tq and 7 at a redshift 2 = 3. The values 
of 7 are indicated in the panel. The error bars assume 30 per cent 
uncertainty on (A) per 5000 km s~^. The line is a least squares 
fit to the points, demonstrating the scaling T1.2 oc {A)~^. 



Cp(wi) 



E 



C{v 



1,1^2 



(5) 



For a given velocity, Cp{v) is a measure of the fraction of 
spectrum that has a significantly different temperature. The 
level of significance of a given value Cp{v) is obtained by 
performing the same analysis on all the randomised spec- 
tra, from which we can compute the probability distribu- 
tion P{Cp). This concludes the description of the statistical 
analysis. 
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Figure 11. Statistical significance that the temperature in a win- 
dow of size 5000 km s"^ starting at velocity v differs from the 
spectrum as a whole. Outlined histogram (left hand scale) is sig- 
nificance in per cent for cold regions, filled histogram (right hand 
scale) for hot regions. Dotted line (left hand scale) denotes the 
flux, smoothed over 5000 km s" Both the hot and the cold re- 
gion in SI (top panel) are assigned a significance of 100 per cent. 
The cold gap in spectrum S2 has a significance of 95 per cent. 
Finally, no significant region is detected in the single tempera- 
ture spectrum S3 (bottom panel). The bottom panel shows the 
A statistic for S3, which can be compared with the corresponding 
curve for SI in Fig. p 



4.3 Temperature calibration 

The method described so far identifies regions of signifi- 
cantly different temperatures. We can get an idea of how 
much the temperatures differ by examining in more de- 
tail the probability distribution P{A) of wavelet amplitudes. 
The shape of this distribution depends mainly on its average 
{A), as pointed out by Zaldarriaga (2001), and demonstrated 
in Fig. y. There is an additional dependence on the slope 7 
of the T — p relation, in that models with larger 7 have a 
higher fraction of wavelets around the mean. 

In turn, {A) correlates strongly with the temperature 
Ts = r5^~^ of the gas, at some overdensity 5 = p/{p)- 
For the mean absorption appropriate for 2: = 3, 5 ~ 1.2. 
Figure ^ shows that Ts oc {A)~^ , for simulated models with 
a range of values for To and 7. So without calibrating the 
{A) — T relation with simulations, one can still estimate 
temperature ratios. If one is willing to normalise the relation 
using simulations, then it is possible to measure Ti.2 directly. 

In the linear regime, the baryonic overdensity 5igm is 
smoothed with respect to the dark matter field 5dm accord- 
ing to SiGM{k) oc S-DM{k)/k^T, at high wavenumbers k (e.g. 
Bi and Davidsen 1997). The wavelet amplitude {A) is a mea- 
sure of the r.m.s. of 5igm on small scales, which explains the 
scafing (A) oc (^i^gm)'^' oc T"!. 



The mean (A) can be directly computed for a given 
spectrum, and the linear relation of Fig. hd can then be used 
to estimate T1.2. In the simulations we measured a variation 
of ~ 30 per cent in (^4) over regions of dv — 5000 km s~ , 
in single temperature models. The corresponding error bars 
are indicated in Fig. |l^. When applying this linear scaling 
to data of velocity extent dv, we will also assume a similar 
error ~ 30(9'i;/5000kms~^)~^''^ per cent. The values 7=1 
and 7 = 5/3 are likely to span the range encountered when 
applying this calibration to real data. 

4-3.1 Summary and illustration 

The method we have devised consists of three steps. 

(1) Compute the contribution of a narrow wavelet to the 
derivative of the spectrum. We used simulation to show that 
the amplitude A of such a wavelet anti-correlates strongly 
with the temperature of absorbing gas (Fig. |l|). The cumu- 
lative distribution P(< A) can be used to characterise the 
p — T relation. For example, we used the maximum distance 
A between the cumulative distributions of two stretches of 
spectrum to judge which stretch corresponds to higher tem- 
peratures (Fig. H) . 

(2) Generate random spectra from the original data by 
scrambling the line list. These random spectra are used to 
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Figure 9. Dependence of the probability distribution P{A) of wavelet amplitudes on To and 7. Panel (a): P{A/{A\)) versus A/{A\) for 
models with logj^gTo = [4.0,4.1,4.2,4.3] (full, dotted, short-dashed and long-dashed lines respectively) and 7 = 1.4. The amplitudes have 
been scaled by the mean wavelet amplitude (Ai) of the coldest model. In panel (b), the distribution P{A/{A)) is scaled using the mean 
of the distribution. In panel (c) we show the normalised distributions P(A/(A)) for models with (log^QTo,7) = (4.18, 1) and (4.18,5/3) 
(lower curves, full and dashed line respectively, off-set vertically by -0.1) and (log]^o^0,7) = (4.34,1) and (4.34,5/3) (full and dashed 
lines respectively, off-set vertically by 0.1). Panel (a) shows that P{A) depends on To, with hotter models having a smaller fraction of 
large wavelet amplitudes, as expected. Panel (b) demonstrates that the difference in shape is mostly due to the difference in mean {A) 
of the distribution. Finally, panel (c) shows that 7 has an effect on the shape of P{A/ {A)) as well, such that models with a larger 7 are 
larger around the mean of the distribution. 



compute the probability distribution P(A) for spectra with- 
out temperature fluctuations. 

(3) Identify unusual regions in the spectrum, that corre- 
spond to large values of |Q|, where Q{vi,V2,dv) is defined 
in Eq. (^). Use a clustering analysis to better characterise 
the size of these regions. Project the resulting clusters on 
to the spectrum to obtain Cp{v) (Eq. ^) and use the ran- 
domised spectra to associate the statistical signiflcance. 

An estimate of the temperature ratio of different regions 
can be obtained by comparing the mean wavelet amplitudes 
and using the relation T oc {A)^^ . 

The final result of applying the above procedure to the 
simulated spectra SI, S2 and S3, is shown in Fig. For each 
value of Cp{v), we flnd the fraction / of random spectra 
that have a comparable cluster, and plot the signiflcance 
100(1. — /) in per cent (or 100(/ — 1.) for cold regions), 
both for hot regions (filled histogram, right hand scale) and 
cold regions (histogram, left hand scale). The hot and cold 
regions in spectrum SI both have a signiflcance of ~ 100 per 
cent, meaning that none of the (200 in this case) randomised 
spectra contain clusters that large. The cold gap in spectrum 
S2 is signiflcant at the 95 per cent level. 

The different p — T relations in the spectra SI and S2 
were imposed on the simulations in post processing. This is 
likely to underestimate the differences in line widths, com- 
pared to models in which the temperature is different during 
the simulation as well. For example, Theuns et al. (2000) 
showed that line widths for simulation with identical im- 
posed p — T relations are measurably different if the un- 
derlying simulations had different temperatures. The reason 
is that simulations with different temperatures vary in the 
amount of Jeans smoothing and pressure induced peculiar 
velocities, which also contribute to the line widths. These ef- 



fects are not captured by changing the p — T relation, which 
only inffuences the thermal broadening. 

Finally, spectrum S3 drawn from a single temperature 
model does not contain any regions signiflcant to more than 
10 per cent. Recall that this simulated spectrum has a large 
jump in optical depth, with mean fluxes F ~ 0.64 and ~ 0.60 
for the first, respectively second half of the spectrum. This 
latter test demonstrates that the procedure of randomising 
spectra works well — the original spectrum is equivalent to 
the randomised spectra for the statistic we investigate here. 
Note also that the method is not very sensitive to the mean 
effective optical depth. We apply the analysis procedure to 
data in the next Section. 



5 THERMAL STATE OF THE IGM 
5.1 Redshift range z = [3,3.6] 

Two of our QSO spectra span the region around z ~ 3.2 
where Schaye et al. (2000) claimed an increase in To. We 
have prepended a stretch of spectrum from QSO 0302-003 
and of 0014+813 to that of APM 0827-1-5255, to have two 
more spectra that straddle that redshift z ~ 3.2. We use 
these spectra to investigate whether our new method inde- 
pendently suggests a change in To around redshift z ^ 3.2. 
For each of the spectra, we have generated 200 randomised 
spectra to judge the significance of the temperature fluctu- 
ations as described earlier. 

QSO 0055-269 does indeed appear to show a very sud- 
den increase in To around z - 3.3 (Fig. |l|). Its low redshift 
half is unusually hot at the 99.5 per cent level, and the high 
redshift half is unusually cold at the 97 per cent level. It is 
worth noting that the low redshift half contain two voids at 



© 0000 RAS, MNRAS 000, 000-000 



Temperature fluctuations in the IGM 11 
QSO 0055-269 QSO 0302-003 and 0827+5255 

— I 1 1 1 1 1 1 1 1 1 1 1 — □ r — I 1 1 1 1 1 1 1 1 1 1 1 1 3 




3 3.2 3.4 3.6 3 3.2 3.4 3.6 

redshift redshift 



Figure 12. Wavelet analysis of the spectrum of QSO 0055-269. 
Top panel: mean flux averaged over a window of dv = 8000 km 
s~^. Bottom panel: temperature measure A{v,dv) (dot-dashed 
line, right hand scale)|.^and statistical significance P{Cp), where 
Cp is defined in Eq. (H), of temperature fluctuations in per cent 
(histogram, left hand scale). The low redshift half of the spectrum 
is unusually hot at the 99.5 per cent level, and the high redshift 
half of the spectrum is unusually cold at the 97 per cent level. 
The jump in temperature appears to be very sudden and occurs 
at a redshift z ~ 3.3. The mean absorption F does not appear to 
undergo a similar strong evolution. 



0.8 



-100 t-i L 



QSO 1422+231 



T 



T 




3.2 3.4 
redshift 



3.6 



Figure 13. Same as Fig.|l|for QSO 1422+231, but for a window 
size of dv = 5000 km s~^. Below z 3.05, the spectrum is unusu- 
ally hot at the 97.5 per cent level. Above z ^ 3.4 it is unusually 
cold at the 99.5 per cent level. In between, there appear to be 
large fluctuations in A(d, dv), which however are not statistically 
significant, when compared to the spectrum as a whole. 

z ~ 3.11 and 3.28 of sizes ~ 17 and 15 Mpc respec- 
tively (deceleration parameter qo = 0.1, Kim et al 2001b). 
This sudden transition looks very similar to what we had for 
our simulated spectrum SI (but there the low redshift half 
of SI was cold). The identification of a statistically signifi- 
cant change in To does not require the use of simulations. 
However, once we need to decide by how much To changes, 
simulations are needed in order to calibrate the relation be- 
tween the wavelet distribution and (To, 7). 

QSO 1422-1-231 also shows an increase in To, albeit not 
so suddenly (Fig. As for QSO 0055=269, the top (bot- 
tom) half is unusually cold (hot), at the 99.5 (97.5) per cent 
level. Note that there are large fluctuations in A in between 
the hot and cold regions, in the interval z — [3.1,3.4]. In 
particular, it would be tempting to identify the hot region 
around z ~ 3.37 as being heated early, and/or the cold re- 



Figure 14. Same as Fig. |l3| for QSO 0302-003 (interval z = 
[3,3.27] and QSO APM 0827+5255 (interval z = [3.27,3.7]). The 
spectrum below z ~ 3.35 is hot at the 100 per cent level, and 
above z ~ 3.35 it is cold at the 99.5 per cent level. 
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Figure 15. Same as Fig. but for QSO 0014+813 (interval 
z = [3,3.2] and QSO APM 0827+5255 (interval z = [3.2,3.7]). 
The spectrum below z ~ 3.15 is hot at the 100 per cent level, and 
above z ~ 3.35 it is cold at the 99.5 per cent level. 



gion around z ~ 3.2 as remaining cold slightly longer. How- 
ever, the regions around z ~ 3 and z ~ 3.5 are much larger, 
hence far more significant. When randomised spectra are 
generated from this spectrum, many will actually contain 
large cold/hot regions as well, making the apparent fluctu- 
ations around z ~ 3.2 not statistically significant. 

The combined spectrum of QSOs 0302-003 (from z = 
[3,3.27]) and APM 0827+5255 (from z = [3.27,3.7]) is 
shown in Fig. |l^. There is a sudden change in tempera- 
ture at 2 = 3.3, with a significance of 100 per cent for the 
low redshift hot half, and 99.5 per cent for the high redshift 
colder region. QSO 0302-003 contains a well known void 
at z ~ 3.17 (Dobrzycki & Bechtold 1991) which is visible 
in Fig. ^ as a region of low absorption. We note that the 
higher resolution data of Kim et al (2001b) suggest a smaller 
size for the void than the earlier data. In any case, compar- 
ing the mean absorption with the temperature measure A, 
it is clear that the high inferred temperatures do not just 
result from the presence of these voids, although they might 
contribute to it. 

The combined spectrum of QSOs 0014-1-813 (from z — 
[3,3.2]) and APM 0827-^5255 (from z = [3.2,3.7]) is shown 
in Fig. |l^. Below z = 3.15, this spectrum is hot at the 100 
per cent level, above z — 3.35 it is cold at the 99.5 per cent 
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Figure 16. Probability distribution of wavelet amplitudes, 
P(j4), for the combined spectrum of QSOs 0302-003 and 
APM 0827+5255, and for QSOs 1422+231 and 0055-269. 
The high redshift halves of the spectra {[3.37, 3. 6], [3. 45, 3.5] 
and [3.35, 3.5] respectively) and the low redshift halves 
([3.05, 3. 3], [2. 95, 3. 05], [3.05, 3.25]) are shown separately. P{A) for 
spectra at the same redshift are very similar, but there is a differ- 
ence between the probability distributions at different rcdshifts. 



level. In between, there are large scale fluctuations similar 
to what we obtained for QSO 1422+231. 

The probability of wavelet amplitudes, P{A), is plot- 
ted in Fig. |l^ for the unusual regions identified in these 
spectra. Consistent with the above evidence of a temper- 
ature jump, we find a change in shape of P{A) around 
z ^ 3.3. Comparing the normalised distribution, P{A/ {A)), 
with the simulations, we find good agreement between the 
high redshift halves and the simulation with (logTo,7) = 
(4.18,5/3), and the low redshift half and the isothermal 
model (logro,7) = (4.34,1) respectively (Fig. |l^. Using 
these values for 7, we can apply the calibration between {A) 
and T1.2 = Tol.2^~^ from Fig.|lC| to obtain an estimate for 
To. The result is in excellent agreement with the tempera- 
ture evolution determined by Schaye et al. (2000) using the 
h~ N cut-off (Fig. Ili). Above redshift z ~ 3.4, the temper- 



ature To ~ io''-^±"-^=', and below that To ~ lo^-SiO-is, an 
increase of around 60 per cent. 

We conclude that all four spectra show evidence for 
a change in temperature, significant at the ~ 99 per cent 
level, over a relatively narrow redshift interval around z ~ 
3.3±0.2. The temperature increase is 60 ± 14 per cent. Such 
an evolution is clearly not consistent with photo-heating of 
a highly ionised, expanding medium, in which case the tem- 
perature would smoothly decrease towards lower z, but is 
consistent with a sudden entropy injection in the IGM re- 
sulting from Hell reionization. 



5.2 Redshift range z = [2.6, 3.2] 

We have combined the four intermediate redshift QSOs 
0636+680, 0956+122, 0302-003 and 0014+813 into one 
longer spectrum analysed in Fig. The lowest redshift 
part of 0302-003 has a cold region of moderate significance 
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Figure 17. Same as Fig. |l^, but for the distribution of wavelet 
amplitudes in units of the mean amplitude. A/ {A)P{A/ [A)). 
Dashed lines refer to the high redshift parts of the spectrum, full 
lines to the low redshift halves. Filled triangles and filled squares 
refer to simulations with {logT(),j) = (4.18, 5/3) and (4.34, 1) re- 
spectively. The high redshift halves are well represented by the 
colder model with a steep T — p relation, the lower redshift halves 
are better fit by the hotter, isothermal model. 



(P ^ 85 per cent). In addition there is a hot region of size 
~ 10* km s~^ at a redshift z ~ 3.1 in the spectrum of 
0014+813. This hot region is also significant if we analyse 
the spectrum of 0014+813 on its own (Fig. 

The wavelet amplitudes and the spectrum itself are 
compared for this peculiar region in the spectrum of QSO 
0014+813 and two other stretches at comparable redshifts, 
in Fig. What makes the spectrum in the top panel pecu- 
liar from the other two, is that lines of comparable strength 
are broader and consequently produce smaller wavelet am- 
plitudes. The analysis with 200 randomised spectra puts the 
statistical significance of the hot region at the P = 99.5 per 
cent level. Using the mean wavelet amplitude to estimate the 
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Figure 21. Flux (full lines, left hand scale) and wavelet amplitudes (filled histograms, right hand scales) for three regions of length 10* 
km taken from the spectrum shown in Fig. The top panel corresponds to the hot region around z ~ 3.1 (v ~ 1.1 X 10^ km s~^), 
the middle and bottom panel correspond to stretches of spectrum at comparable redshifts, that are not deemed unusual. Note the small 
fraction of narrow lines in the hot region, as compared to the other two. 



temperature, we find that the difference in To is 60 ± 30 per 
cent. This higher temperature is presumably a consequence 
of the region undergoing reionization later than the rest of 
the spectrum, although other processes might contribute as 
well. 



5.3 Redshift range z = [1.5, 2.4] 

We have plotted the combined spectrum of the four low 
redshift QSOs 0515-4414 (z = [1.5,1.7])], 2233-606 {z = 
[1.8,2.2]), 1122-1648 {z = [1.9,2.4]) and 2217-2818 (2 = 
[1.9,2.3]) in Fig. ^ QSO 2217-2818 contains a region of 
size ~ lO** km s~ 100 Mpc comoving size in a 

{Q,m,^A) = (0.3,0.7) cosmology) which is colder at the 99 
per cent significance level, when compared to 200 random 
realizations. The mean absorption does not appear to be 
unusual in this region. When analysing the spectrum of QSO 
2217-2818 on its own, the cold region is equally significant 
at the 99.5 per cent level (Fig. ||). 



The wavelet amplitudes and the spectrum itself are 
compared for this cold region and two other stretches at 
comparable redshifts, in Fig. |2^. The differences between the 
top panel and the two others, is similar to what we found for 
the simulated spectrum SI in Fig. 0, with the narrower lines 
in the top panel generating larger wavelet amplitudes. Using 
the mean wavelet amplitude to measure the temperature, we 
find the difference in To is 77 ± 36 per cent. 

At these lower redshifts, the number of metal lines, and 
the fraction of lines blended with metal lines, is no longer 
negligible. In order to minimise the possibility that the de- 
tected cluster of narrow lines is significantly contaminated 
by such lines, two of us have independently fit the spectrum 
of QSO 2217-2818 using VPFIT . For each line not directly 
identified as a metal line, we have checked whether there 
is a consistent detection of Lyman- /3 as well (the data are 
noisier at Lyman-7, but we also demand consistency at that 
transition where possible, i.e. for z > 2.137). Unfortunately, 
metal transitions often have multiple components, and so if 
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Figure 24. Flux (full lines, left hand scale) and wavelet amplitudes (filled histograms, right hand scales) for three regions of length 10* 
km s~^ taken from the spectrum shown in Fig. The top panel corresponds to the cold region around z ~ 2.1 (v ~ 1.5 X 10^ km s~^), 
the middle and bottom panel correspond to stretches of spectrum at comparable redshifts, that are not deemed unusual. Note the large 
fraction of narrow lines and the associated high wavelet amplitudes in the top panel. The two other panels do not contain such narrow 
lines. 



a transition were to be missed, it can contribute significantly 
to the signal. We have performed the wavelet analysis with 
these two independent line lists, and the results are nearly 
identical. Both detect an unusually large wavelet signal in 
the same region, at significance levels of 99.5 and 98.5 per 
cent respectively. 

From a theoretical point of view, if the temperature 
difference were due solely to a difference in reionization red- 
shift, then a region could have an unusually low temperature 
if it had been reionised very early, or not at all. Assuming 
that T oc (1 + zY'^ after reionization, the early redshift 
would have to be ^ 4.5. Both these possibilities seem un- 
likely to us. Given that we only detect one such an unusual 
region over this redshift range, it seems more plausible that 
a combination of several effects, such as residual and par- 
tial contamination by metal lines, the presence of large scale 
structure, and the gradual temperature decrease due to the 
expansion of the universe, all contribute to making this re- 



gion appear unusual. In addition, the voids in the spectrum 
at redshifts z ~ 1.9, 2.2 and 2.3 (Kim et al 2001a) will con- 
tribute to the signal as well. 



6 DISCUSSION 

Schaye et al. (2000) determined the evolution of the p ~ T 
relation of the IGM by measuring the cut-off in the b — N-^^ 
distribution of 9 high resolution QSOs. (Only two of those 
are in common with the ones used in the present analy- 
sis). They used high resolution hydrodynamical simulations 
to calibrate the relation between the b — A'jjj cut-off and 
the underlying p — T relation (Schaye et al. 1999). They in- 
ferred relatively high temperatures. To ~ 10^'^K, 7 ~ 1.3 
at high redshifts z = [4.5,3.5], an apparently sharp rise to 
To ~ 10'* ''K around 2; ~ 3 with an associated dip in 7 1, 
followed by a gradual decrease in To and increase in 7 to- 
wards 2 ~ 2. This temperature evolution is not consistent 
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Figure 18. Temperature at the mean density, To, versus redshift. 
Open symbols with error bars are taken from Schaye et al. (2000). 
Filled symbols with error bars and thicker lines refer to estimates 
of To based on the mean wavelet amplitudes and the calibration 
of Fig. |lo| with errors on Tq of 30 per cent for a stretch of spec- 
trum of 5000 km s~^. Errors in the redshift direction refer to 
the redshift extent of the sampled region. Symbol types refer to 
the QSOs (note that APM 0827+5255 and 1422+231 were used 
in the analysis of Schaye et al. 2000 as well.) We have assumed 
7 = 1 below z = 3.3, and 7 = 5/3 above z = 3.3 to convert T1.2 
to Tq. Dotted line is the temperature evolution in a simulation 
where Hen reionizes at redshift ~ 3.4. It fits the data very well. 



with photo-heating in the optically thin limit, as determined 
from the Haardt & Madau (1996) tracks. 

Ricotti et al. (2000) wrote down a partly theoretically 
motivated parametric function that describes the 2D distri- 
bution of lines in the b — A^jjj plane. They employed dark 
matter simulations with an added 'pressure term' to mimic 
the effects of thermal smoothing, to calibrate the fitting pa- 
rameters of the function, in terms of the simulation param- 
eters. Their mock spectra were analysed with AUTOVP 
(Dave et al. 1997) . Applying their method to published line 
lists, they obtained a thermal evolution consistent with that 
deduced by Schaye et al. (2000), although their error bars 
are large. In particular, they also stressed the small values 
of 7 ~ 1 around z ~ 3. 

Bryan & Machacek (2000) used high resolution Eule- 
rian simulations to measure To from the b — A'jjj cut-off. 
These authors stressed that the position of the cut-off also 
depends sensitively on the amplitude erg of the underlying 
dark matter power spectrum. Bryan & Machacek used the 
Voigt profile fitter described in Zhang et al. (1997), which 
fits Gaussians to the optical depth distribution of absorption 
lines without attempting to deblend the absorption lines in 
terms of Voigt profiles, in contrast to VPFIT or AUTOVP . 
Theuns et al. (2000) showed that these other fitters appear 
much less sensitive to as , which might partly explain the dif- 
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Figure 19. Combined spectrum of four z ~ 2.8 QSOs as func- 
tion of velocity v. Top panel; redshift range for each QSO; mid- 
dle panel: mean flux averaged over dv = 5 X 10^ km s~^; bot- 
tom panel: temperature indicator A{v, dv) (dot-dashed line, right 
hand scale) and significance P in per cent (histogram, left hand 
scale). There is a significantly hotter region (P 99 per cent) 
of size ~ lO-* km s-i at ~ 1.1 X 10^ km s"! {z ~ 3.0) in the 
spectrum of QSO 0014+813. 
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Figure 20. Same as Fig. ^ but for the spectrum of QSO 
0014+813, using a window size dv = 5000 km s-l. The hot region 
around z ~ 3 has a significance of 98.5 per cent. 



ference in erg dependence. (Theuns & Zaroubi (2000) demon- 
strated that the wavelet analysis as described here is also 
not very sensitive to ag) Bryan & Machacek inferred high 
temperatures as well, and also concluded that the heating 
by Hell might be underestimated in the traditional models 
based on the optically thin limit. 

McDonald et al. (2000) used a slightly different proce- 
dure to measure To from the widths of the lines. They only 
use lines that are well-fit by a single, unblended Voigt pro- 
file. Tq is determined from the widths of those fits, correcting 
for other broadening contributions as a function of column 
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Figure 22. Combined spectrum of four S ~ 2 QSOs as function of 
velocity v. Top panel: redshift range for each QSO; middle panel: 
mean flux averaged over 9ti = 5 X 10"^ km s~^; bottom panel: 
temperature indicator A(d, dv) (dot-dashed line, right hand scale) 
and significance P in per cent (histogram, left hand scale). There 
is a significantly colder region (P ~ 99 per cent) of size ~ 10* 
km s~^ at I) ~ 1.5 X 10^ km s^^ (z ^ 2) in the spectrum of QSO 
2217-2818. 
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Figure 23. Same as Fig. ^ but for the spectrum of the low 
redshift QSO 2217-2818, using a window size dv = 5000 km s'^. 
The cold region in the middle of the spectrum has a significance 
of 99.5 per cent. 



density, using simulations. They find this correction to be 
large for low density gas, and hence resort to measuring Ts 
at an overdensity of 5 ~ 1.5. They find no evidence for evo- 
lution in Ts - lO'^-^K with redshift, from z ~ 3.9 to 2.4. 
The redshift evolution of the temperature at a given over- 
density, Ts, depends on S, because the photo- heating rate 
depends on density (Miralda-Escude & Rees 1994; see also 
Fig. (3) in Schaye et al. 2000). It is not clear whether the 
measurement of T1.5 by McDonald et al. (2000) is in contra- 
diction with the determination of To by the other groups, 



given the large uncertainty in 7 quoted by all groups. The 
fact that their analysis is based on a single hydrodynamical 
simulation makes it difficult to judge the robustness of their 
analysis procedure. In addition, these authors averaged their 
results over wide redshift bins, which hampers their ability 
to detect sudden changes in To or 7. 

Zaidarriaga, Hui & Tegmark (2000) used the small scale 
cut-off in the power spectrum as a measure of Tq. They find 
similar temperatures as McDonald et al. (2000), but unfor- 
tunately only use dark matter simulations to calibrate their 
results. Also RoUinde et al. (2001) use an inversion proce- 
^dure to measure To which is tested only on low resolution 
dark matter simulations. Such simulations are inadequate to 
study the thermal properties of the IGM. Even simulations 
that include hydrodynamics need to be very accurate in or- 
der to resolve all the line broadening mechanisms properly 
(Theuns et al 1998; Bryan et al. 1999). Therefore this type of 
analysis only shows that the method can deduce the temper- 
ature of the gas painted on top of a dark matter simulation 
using an untested prescription. 

All determinations so far lead to higher values of To 
than were expected for photo-heating by a UV-background 
dominated by QSOs, such as computed by Haardt & Madau 
(1996). One way to increase To for a given UV-background 
is to increase the baryon fraction uii, = Q,bh? , since the heat- 
ing rate is proportional to the physical density of the gas 
(Miralda-Escude & Rees 1994). Theuns et al. (1999) demon- 
strated that higher lj^ makes lines fitted by VPFIT consid- 
erably broader as Tq increases. A higher value of loi, would 
also allow galaxies to make a significant contribution to the 
UV-background, as seems to be implied by the measure- 
ments of Steidel et al. (2001) of the large escape fraction 
of UV-photon from Lyman-break galaxies (Haehnelt et al. 
2001). Without such a high uji,, it would be difficult to ex- 
plain the observed high opacity of the IGM. However, such 
high values of ujt would appear to violate the tight con- 
straints on ujh ~ 0.0193 ± 0.0014 determined from the deu- 
terium abundance (Buries & Tytler 1998; Kirkman et al. 
2000). 

Alternatively, higher values of To can result from other 
heating mechanisms. The currently most popular one is that 
the calculations based on assuming an optically thin gas un- 
derestimate the photo-heating rate, by a factor of up to a 
few (Miralda-Escude & Rees 1994; Abel & Haehnelt 2000). 
Therefore, delayed Hell reionization could be responsible for 
the increase in To around z ~ 3, as suggested by Schaye 
et al. (2000), Ricotti et al. (2000) and Bryan & Machacek 
(2000). In fact, even the constant To favoured by McDon- 
ald et al. (2001) requires an extra heating source, such as 
Hell reionization. 

If the sources responsible for reionizing Hell are bright 
but scarce, then the temperature distribution could be- 
come quite inhomogeneous, with large Hem regions sepa- 
rated from cooler, Hell zones (Miralda-Escude & Rees 1994). 
The ionising radiation would quickly over run the low den- 
sity voids, impulsively heating the gas, while the more over 
dense filamentary regions would remain neutral for longer 
(Miralda-Escude, Haehnelt & Rees 2000). Since the lower 
density gas then will be hotter, its neutral fraction would 
become even lower because of the T~^''' temperature de- 
pendence of the recombination coefficient. As a result, the 
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contrast in absorption between voids and filaments would 
increase even further. 

The large jump in To of 60 ± 14 per cent over the red- 
shift interval z = [3.5,3.1] discussed in Section 5.1 strongly 
suggests a sudden entropy injection, likely resulting from 
Hell reionization. Most of the wavelet signal comes from gas 
at modest over densities 1 ^ 5 3, which produces lines of 
column density ~ 10^^ cm~^. This is because weaker lines 
do not contribute significantly to the wavelets, and stronger 
lines are both fewer, and partly explicitly ignored in the 
analysis. This gas traces the modestly over dense filamen- 
tary pattern seen in simulations, and our analysis suggests 
it becomes reionized in Hell in a fraction ~ 0.2 of a Hubble 
time. This also naturally explains the relatively high values 
of To found at lower redshifts. 

A sudden increase in To influences the gas distribution 
in shallow potential wells. Theuns et al (2000) explicitly 
showed how increased heating can drive gas out of the cen- 
ters of filaments into the lower density surroundings. The 
velocity of this gas also contributes to the line widths as 
determined by Voigt-profile fitting. Such impulsive heating 
may infiuence the formation of halos of virial temperatures 
~ 2. X 10* K (Gnedin 2000b). 

After reionization, the IGM temperature will fall at a 
rate To « (1 -I- z)^'^, slower than pure adiabatic expansion 
because of the continuing photo-heating. Over a redshift ex- 
tent Sz ~ 0.5 of the Lya forest of a QSO, this corresponds to 
a drop of ~ 27 and 32 per cent at redshifts 3 and 2.5, respec- 
tively, between the high and low redshift parts of the QSOs 
spectrum. This will contribute to making the top half of 
QSO 0014-1-813 unusually hot, and the bottom half of QSO 
2217-2818 unusually cold. Similarly, if Hell reionization were 
delayed in one region with respect to another, for example 
as a result of QSO clustering, then we might expect to see 
differences in To of order 40 per cent, given our estimate of 
z = 3.3 ± 0.15 for the duration of the reionization epoch. 
Such entropy fluctuations are likely to remain present for 
a long time afterwards, given that the IGM evolves nearly 
isentropically after reionization. 

Why don't we see fluctuations on smaller scales? The 
analysis of spectrum S2 suggests that the present method 
can detect fluctuations of amplitude 50 per cent in To on 
scales as small as 5 x 10"^ km s~^. However, our simulations 
are clearly somewhat idealised: they lack large scale power, 
and in addition have an imposed p — T relation without any 
scatter. In the real universe, the presence of large scale power 
and scatter in the p — T will tend to increase the intrinsic 
fluctuations in the wavelet amplitudes, making it harder to 
recognise a stretch of spectrum with truly different To or 
7. This is particularly true around reionization, where small 
scale fluctuations would be hard to detect in the presence 
of a large jump in the spectrum as a whole. However, the 
absence of small scale fluctuations away from reionization 
suggests that the IGM does obey a reasonably well deflned 
T — p relation on sufficiently large scales > 5000 km s~^. 
Given the arguments above it also suggest that reionization 
occured relatively well synchronised throughout the IGM. 



7 SUMMARY 

We have presented a method that uses wavelets to char- 
acterise objectively the line widths of Lya lines in a QSO 
spectrum. Colder gas produces narrower absorption lines, 
for which the wavelet amplitudes are larger. We described 
a statistical method, based on recognising regions were the 
distribution of wavelet amplitudes is unusual, to identify re- 
gions where the temperature at the mean density. To, dif- 
fers significantly from the mean. The significance level is 
determined by performing the wavelet analysis on spectra 
obtained from the data after scrambling the list of absorp- 
tion lines, thereby destroying all correlations imprinted by 
(temperature) fiuctuations. 

We applied this method to mock spectra generated from 
high resolution hydrodynamical simulations. We imposed 
several p — T relations on the output from these simula- 
tions, and used these to construct spectra with intrinsic tem- 
perature fiuctuations. We demonstrated that the method is 
able to detect fluctuations in To of ~ 50 per cent in regions 
as small as 5000 km s"'^ at redshift z = 3, at the 95 per 
cent level. We also showed that by calibrating the relation 
between the average wavelet amplitude, (A), and To using 
our hydrodynamical simulations, one can estimate the IGM 
temperature with an rms error of order of 30 per cent, per 
stretch of spectrum of 5000 km s~^. 

We have applied the method to eleven high resolu- 
tion QSO spectra. There is strong evidence of a sharp in- 
crease in To by a factor ~ 60± 14 per cent over the red- 
shift interval z = [3.5,3.1]. We interpreted this as evidence 
for Hen reionization, and concluded that it took of order 
5z ^ 0.4 to reionize the filaments at overdensity ~ 3 that 
our method is most sensitive to. Our estimate for the tem- 



peratures above and below the jump are To 



10* 



and To ~ 10* 



% respectively. These values are in excel- 



lent agreement with those obtained by Schaye et al. (2000), 
who used a very different method. 

At lower redshifts, we find that the IGM appears to 
follow a reasonably well defined T — p relation on scales 
larger than 5000 km s"'^. In 8 QSO spectra, we only found 
two unusual regions, significant at the more than 99 per cent 
level. QSO 00144-813 has a hot region of size 10* km s"^ at 
z ~ 3.1, where the temperature is around 60± 30 per cent 
higher. This could be a relic of Hell reionization, enhanced 
by the over all descrease in To along the spectrum as a result 
of adiabatic expansion. QSO 2217-2818 has a cold region of 
similar size, at z ~ 2.3, where To is lower by ~ 77± 36 per 
cent. Again in addition to adiabatic expansion, the three 
voids in this spectrum will contribute to making this region 
unusual. 
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